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Abstract 

7~J We study by means of an Eulerian-Lagrangian model the statistical properties of velocity and acceleration of a neutrally-buoyant 
finite-sized particle in a turbulent flow statistically homogeneous and isotropic. The particle equation of motion, beside added mass 
and steady Stokes drag, keeps into account the unsteady Stokes drag force - known as Basset-Boussinesq history force - and the 
i— I non-Stokesian drag based on Schiller-Naumann parametrization, together with the finite-size Faxen corrections. We focus on the 
case of flow at low Taylor- Reynolds number, Re,\ -31, for which fully resolved numerical data which can be taken as a reference 
£^ are available (Homann & Bee 651 81-91 J. Fluid Mech. (2010)). Remarkably, we show that while drag forces have always minor 
effects on the acceleration statistics, their role is important on the velocity behavior. We propose also that the scaling relations for 
' O the particle velocity variance as a function of its size, which have been first detected in fully resolved simulations, does not originate 
^2 from inertial-scale properties of the background turbulent flow but it is likely to arise from the non-Stokesian component of the drag 
rj~] produced by the wake behind the particle. Furthermore, by means of comparison with fully resolved simulations, we show that 
^ the Faxen correction to the added mass has a dominant role in the particle acceleration statistics even for particle with size in the 
Q inertial range. 
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1. Introduction 

The exact dynamics of a material particle in an inhomoge- 
neous unsteady flow involves nonlinear equations that can be 
treated analytically only in approximate for m 1,2 - 3 ' 4 . For this 
reason, several simplified models for the hydrodynamic forces 
acting on a particle have been proposed in the literature^. It 
is still unclear however to what extent these models provide 
an accurate description in turbulent flow conditions, even in an 
averaged or statistical sense. A proper statistical description of 
particle dynamics would be a first important step towards build- 
ing constitutive equations for the particulate phase carried by 
turbulent fluids. It would also be of practical importance for the 
many environmental phenomena and industrial applications in 
which particle suspensions in turbulence are involved. We in- 
tend here to contribute to this goal by carrying out refined sim- 
ulations of particles in turbulent flow, and discriminate whether 
the particle model employed leads to physically sound results, 
in agreement with recent experiments and with fully resolved 
direct numerical simulations. 

In previous studies, we have addressed the dynamics of small 
material particles with a description based on a minimal La- 
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grangian model accounting for pressure gradient, added mass 
term, and steady Stokes drag force 7,8 . While this system pro- 
duces several features of particle dynamics, like clustering and 
segregation as well as single-/multi-time statistics of accelera- 
tion and velocity, it fails to predict some statistical properties, 
particularly when the size of the particle is progressively in- 
creased above the dissipative scale of turbulence^. In order to 
better understand these discrepancies we have focused on the 
case of finite-sized and neutrally buoyant particles. We have 
proposed that Faxen corrections are the essential ingredients to 
account for the statistical properties of finite-sized particle ac- 
celeration in turbulence— Numerical predictions of the particle 
acceleration variance (and of its probability density function) 
based on the Faxen argument were compared with experiments 
made in wind-tunnel— and resulted in agreement with Von Kar- 
man flow measurements-is. In a more recent work, experimen- 
tal measurements support also other trends highlighted by the 
Faxen model: the effect of decrease of acceleration flatness as 
a function of the particle-size and the corresponding growth of 
the correlation time of the acceleration—. 
Another series of numerical studies were recently conducted by 
Homann and Bee— (HB in the following). These authors em- 
ployed a direct numerical approach. They tracked the motion 
of a neutrally buoyant finite-sized particle in a turbulent flow 
by enforcing the no-slip velocity at the particle surface via a 
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penalty method on the discretized Navier-Stokes equation. In 
such a way they have been able to access the dynamics of a sin- 
gle finite-size particle in the diameter range [2, 16]rj in a moder- 
ately turbulent flow at Re,\ = 32. Both velocity and acceleration 
statistics were investigated. Therefore, HB measurements pro- 
vides a set of reference data against which one can test particle 
Lagrangian models, as attempted here. The scope of this work 
however goes beyond the validation of a model equation. We 
aim also at having a physical picture of the statistical dynam- 
ics of particles. We have specific questions in mind: What is 
the statistical effect of the drag, particularly the trailing wake 
drag, on the dynamics of a neutrally-buoyant finite-size parti- 
cle in turbulent flow? Does it modify the acceleration statistics 
or rather the velocity one? Is the role of Faxen correction still 
relevant for particles with size in the inertial range? We will 
see how our study provides an answer to these questions and a 
possible interpretation of the phenomenological picture. 
The paper is organized as follow. In section 2 we describe the 
approach adopted in this study, we introduce the Lagrangian 
modeling of the particle dynamics, its numerical implementa- 
tion, and some expected trends for the particle velocity and ac- 
celeration in the vanishing-size limit. In section 3 we present 
the results of the numerical study, starting from the particle 
Reynolds number behavior, and addressing then acceleration 
and velocity statistics as a function of the particle size. Com- 
parison with direct numerical simulation data is discussed in 
detail in section 4. In the conclusions we summarize the main 
results and and give suggestions for possible future experimen- 
tal/numerical investigations. 



continuum fields Du/Dt and u over a sphere of radius r p cen- 
tered at the particle position, respectively: 

Z)ul -, i C Du r> 

— = (4/3*7*)-' — (x,t)d 3 x (4) 
L Dt iv Jv Dt 

[uls = (Ant^T 1 §u{x,t)d 2 x (5) 

The history force is based here on the Basset-Boussinesq dif- 
fusive kernel, ~ (t — t) -1 ' 2 , while f/, is the time over which 
the memory effect is significant. The non-Stokesian drag coef- 
ficient CR e models the effect of the drag induced by the pres- 
ence of a wake behind the particle. Of course in a Lagrangian 
model of particle dynamics, which is only one-way coupled to 
the fluid flow, no wake can be produced. Therefore, we resort to 
a model: The well known Schiller-Naumann (SN) parametriza- 
tion— . The CR e coefficient, which is a function of the particle- 
Reynolds number based on the diameter size d p = 2r p and on 
an estimator of the slip velocity, Re p = \[u] s - v| d p /v, is cho- 
sen to have the form CR Cp = 0.15 • Re° p M1 considered to be a 
good approximation whenever Re p < 1000^. We note also 
that direct numerical simulations of the flow around a solid par- 
ticle maintained fixed in a turbulent flow shows a good agree- 
ment between the real force acting on the particle (as com- 
puted from strain tensor at the surface of the particle) and the 
drag computed from the slip velocity with Schiller-Naumann 
parametrizationiiii. 

2.2. Numerical implementation 



2. Methods 

2.1. Particle equation of motion 

We consider a Lagrangian equation of motion one-way cou- 
pled to a continuum flow u = u(x, t). Such an equation 
keeps into account the pressure gradient and added mass term 
(~ Du/Dt), the drag force and the volume and surface Faxen 
corrections. The drag force is divided into three parts: the 
steady Stokes drag, the unsteady Stokes drag force or History 
force, and the non-Stokesian drag force. All together it reads as 
follow: 




where r p is the particle radius, v the kinematic viscosity, f3 the 
density coefficient /3 = 3p//(p/ + 2 p p ). Following— the Faxen 
corrections are expressed as volume and surface average of the 



We aim at studying the statistical signature of the different 
forces acting on the particle. For this reason in our numerical 
simulations we follow the trajectories of four species (or fami- 
lies) of particles with slightly different evolution equations. The 
first family is described by ([TJ, it includes only the added-mass 
term and the steady Stokes drag and their Faxen corrections. It 
will be called Faxen model with Stokes drag. The second fam- 
ily is defined by (Q~|i+(|2]i, hence it includes also the history force. 
The third is based on (Q~|i+([3]l, therefore the Schiller-Naumann 
correction is here included but not the history force. Finally the 
forth family ([TJ+fO+fO keeps into account all the effects. 



Particle type 


Eq. 


Faxen model with Stokes drag 

+ History force 

+ Non-Stokesian drag 

+ Non-Stokesian drag + History force 




CDI2> 
•EH) 



Table 1 : The four particle families considered in the present study and their 
respective equations of motion. 

In total we integrate simultaneously the dynamics of N p = 
1.28 ■ 10 6 particles. This ensemble is divided into 4 fami- 
lies, each family having 8 different sizes, in the diameter range 
d p G [0, 32] 77. This amount to 3.2 ■ 10 4 particles per type. 
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2.2.1. Faxen Forces 

The implementation of Faxen averages in our simulation is 
based on the Gaussian approximation proposed in a previous 
study 10 . The volume average of fluid acceleration at parti- 
cle position is replaced by a local interpolation at the particle 
position of the continuum field after convolution by a three- 
dimensional Gaussian envelope G(x), with unit volume and 
standard deviation <x. Convolutions are efficiently computed 
in spectral space, the volume averaged field hence reads: 



Du 

Dt y 



f) 



r goo - 

Ji? Dt 



(x - x', t) d x' 



G(k)g(M) 



(6) 



(7) 



where TfFT ~ denotes a discrete inverse Fourier transform on 
a cubic grid with A^ 3 nodes, while the over script ~ indicates 
a direct Fourier transform, G(k) = exp(-cr 2 k 2 /2) being the 
Fourier transform of G(x). We note that by setting the stan- 
dard deviation <x = r p j V5, in the limit of small radii one gets 
G(k) a 1 - (r?/ 10)k 2 + 0(r p ) which leads to the correct first or- 
der Faxen correction in real space, i.e., u + (r 2 /10)Au + 0(r p ). 
Analogously the surface average reads: 



[u(XJ)]s = hJrS** [u(x,f)]y ) 



(8) 



= OTT _1 [(l - cr 2 k 2 /3) G(k) u(k, f)] . (9) 

For clarity in Figure Q] the shape of the two convolution ker- 
nels (volume and surface) in real space are shown. The figure 
also shows for comparison the so called optimal convolution 
kernels, corresponding respectively to to a three dimensional 
spherical gate function of volume 4/3 nr p for volume average 
and to a delta function over a spherical shell for surface average. 
Of course the implementation of such optimal averages in real 
space would be computationally more expensive. In section [4] 
we will investigate in detail the bias induced on the particle dy- 
namics by the use of such a Gaussian approximation for Faxen 
corrections instead of the rigorous definition. 

2.2.2. History Force 

The Basset-Boussinesq history force can be computationally 
very expensive. This is due to the fact that the integral which 
is involved should be performed at each time-step on the full 
particle history. Furthermore the diffusive kernel, ~ (t - t) -1 ' 2 
has a very slow decay and require a long memory time - virtu- 
ally f/, = oo - to reach convergence. It is known however that 
the diffusive kernel overestimate the history force for particles 
characterized by finite Reynolds numbers Re p The formation 
of a trailing wake either stationary, unstationary or even turbu- 
lent is always associated to history kernel which decay faster 
than the Basset-Boussinesq—. Equivalently we can say that in 
Re p » 1 conditions the Basset-Boussinesq history force should 
have a shorter memory time windows f/,. This latter idea have 
been exploited in the computational approach called windows 
method—. Recently a more performant method based on the 
fit of the diffusive kernel tail via a series of exponential func- 
tions has also been proposed—. In this study we adopt a simple 
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Figure 1: Real space one-dimensional projection, on the direction 
(x,0,0), of filter functions. The volume gaussian filter is G(x) = 
(l/( ^/2ncr)) 3 exp(-x 2 /(2<x 2 )), while the surface convolution kernel turns out 
to be 5(x) = (x - /(2o" ))G(x). The optimal shape of volume and surface filter 
functions are also shown. The volume filter function is a three dimensional 
spherical gate function with total volume Aj'i,nr i p , while the surface one is a 
delta function over the surface of the sphere but normalized in such a way that 
its volume integral is 4?rr 2 . 



windows approximation: instead of setting tj, = oo we chose 
th — 10t, ; but keep the diffusive kernel functional form. This 
choice is based on the observation that in the turbulent condi- 
tions considered in our study after a time 10r^ the Lagrangian 
signal c/duJs - \)/dt is already completely uncorrelated. The 
short memory on the history force is therefore not given by the 
specific kernel form (which is indeed almost flat in our case) 
but from the relatively short correlation time of the turbulent 
flow. With this choice, f;, corresponds approximately to 10 3 
time-steps of our simulations, which are stored and used for 
the discrete estimation of the history integral at each time-step. 
We note that our f/, satisfy the criterion given in ref.21& for the 
windows method and it has a double extension in time steps 
respect to the time windows considered in the numerical val- 
idations considered in ref.— . We have also performed a pos- 
teriori check in which the pre-recorded af([u] s - \)/dt signal 
has been used to compute the history force with different val- 
ues of the windows length. This test has further confirmed the 
convergence and reliability of the adopted implementation. 

2.2.3. Eulerian dynamics 

A suitable turbulent flow is generated by integrating the 
Navier-Stokes equation in a cubic box of size L — 2n with pe- 
riodic boundary conditions. The flow is forced on the largest 
shells in spectral space, on the wave-vectors for which the con- 
dition k 2 < 2 2 (2ttL) 2 is satisfied. The force we adopt in this 
study keeps fixed the amplitude of kinetic energy of the large 
scales. More details concerning the values of relevant input and 
output quantities of the numerical simulation of this turbulent 
flow are provided in Table [2] 

2.3. Faxen corrections and small particle limit predictions 

In the study of Homann and Beoii a derivation of the func- 
tional behavior of the variance of particle velocity and accel- 
eration in the limit of vanishing particle diameters d p has been 
proposed. The argument is based on a perturbative expansion 
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of the Faxen correction for the velocity. This reads asQ: 



[u] s -u+^Au + 0(4). 



(10) 



Furthermore, the hypothesis of a spatially homogeneous par- 
ticle distribution in the limit d p — > is made. By squaring 
Eqn. dTOb . retaining only quadratic terms in d p , and averaging 
over the particle ensemble and in time, (. . .), one gets: 

^>=.<«v|<«A»> = («v|2 = <«v|(|) ! ,ai) 

where s = (v/2)<(Vu + (Vu) r ) 2 ) = v Zr 3 JT 3 uAu d 3 x is the 
mean energy dissipation rate and A = (5v(u 2 )/e)^ 2 is the Tay- 
lor micro-scale. If we instead (i) derive with respect to time 
Eqn. (fTOb . (ii) make the assumption D/Dt d/dt, and then (iii) 
square and average the result, we obtain an approximate predic- 
tion for the acceleration variance 

dl 



D(Vu) 



Dt 



(12) 



We will see in the following to which extent these approxima- 
tions can be considered appropriate to describe the particle be- 
havior. We note that in the simulations we have direct access 
to the values (([u] s ) 2 ) and (([Du/Dt] v ) 2 ) which can be used 
for comparison. Finally, it is also worth noting that the particle 
Reynolds number Re p is proportional to |v - [u] s |, this means 
that in the small-particle limit the leading order is O(d^), hence 
one expects Re p ~ dl. 









1 


















. „2 J/2 , . 
«U >> dp/v 








































/ jr 




















/ p'/ 




















/ / ..- A 












Faxen Model with Stokes drag 










+ Non-Stokesian drag 


— 






+ History 










+ Non-Stokesian drag + History 










fitted slope 4.50 +/- 0.05 





40 



Figure 2: Mean particle Reynolds number, (Re p ), versus the diameter in Kol- 
mogorov scale units, d p /r], form the four different models. A power-law fit to 
the Faxen model with Stokes drag on the smallest particle sizes is shown, we 
get a 4.5 slope. For the same model, the maximal Reynolds number u rms d p /v 
is reached by the largest particles. 



3. Results 

3.1. Particle Reynolds number 

We begin investigating the particle Reynolds number Re p as 
a function of the particle diameter, with measurements reported 
in Fig. [2] First we note that in the range d p € [3.2, 32]/; the 
mean Reynolds number varies considerably, three order of mag- 
nitudes from 10 -1 to about 10 2 . We can immediately observe 
that the particle models which are making use only of Stokes 
drag forces - that is to say based on the assumption Re p < 1 - 
can not be considered entirely consistent. Such models under- 
estimate the actual drag on the particle. This is clearly notice- 
able for the Faxen model with Stokes drag in the large-a^ range, 
when Re p attains the maximal value u rms d p /v, corresponding to 
a ballistic particle velocity v not varying in time and not corre- 
lated to the local fluid velocity [u] s . We note instead that for 
all the models in the small particle limit we have a steeper scal- 
ing (slope 4.50 ± 0.05), close to the expected d p . Hence, in the 
small particle regime v and [u] s are highly correlated, differ- 
ing only by 0(d p ) terms. We also notice that, while the history 
force produces just a shift, the non-Stokesian drag term changes 
the slope in the large particle regime. This apparently minimal 
variations have, as we will see later on, important consequences 
on the statistics of the particle velocity variance. 



1 Note that (10) corrects a typo contained in ref on the numerical coeffi- 
cient in front of rf^Au, which had affected all the predictions proposed in that 
study. Our calculations for eq. lilt provides a coefficient 5/3, instead of 1/100 
given ini4. Furthermore in eq. H2\ we find the coefficient 1/12, and not 1/20. 
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Figure 3: Normalized single-component particle acceleration variance for par- 
ticles of different sizes: <(fli) 2 >/<(D,i(;) 2 > vs. dp/17. Here <(D,«/) 2 > is the fluid 
tracer acceleration variance or equivalently the Eulerian acceleration averaged 
over time and space. The behavior of the four different models adopted is re- 
ported. The dashed-dotted line represent the behavior of particle acceleration 
variance expected to originate from Faxen corrections in the small particle limit 
(see eq 112) ). The dotted line represents the ratio of the variance of the volume 
filtered eulerian field [Du/Dt] v to the fluid acceleration variance. 
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Table 2: Parameters of the numerical simulation: N number of grid points per spatial direction; 8x = 2x/N and 6t are the spatial and temporal discretization; v is the 
value of kinematic viscosity; e the mean value of the energy dissipation rate, rj = (v 3 /e)'' 4 and r,, = (v/e) 1 ' 2 are the Kolmogorov dissipative spatial and temporal 

scales, Urms = ((n/H,-)y/3) ' the single-component root-mean-square velocity, A = (15 v H 2 mJ /e)'' 2 the Taylor micro-scale, 7"£ = (3/2)« 2 mJ /s and = !i,- mI 7£ are 
the Eulerian large-eddy-turnover temporal and spatial scales; = u rms A/ v the Taylor scale based Reynolds number. t m is the total simulation time in statistically 
stationary conditions and the total time-span of particle trajectories. 



3.2. Acceleration statistics 

We examine now some statistical properties of the accelera- 
tion. In Figure[3]we show the behavior of the single-component 
particle acceleration variance (a 2 ) normalized by the fluid ac- 
celeration variance ((D,m,) 2 ) as a function of the particle size in 
dpi r\ units. It is remarkable to note that all the particle models 
leads to very similar results. The History force or the Non- 
Stokesian drag have no effect, at least for this observable. The 
overall trend of the acceleration variance is dominated by the 
Faxen Volume correction, {a 1 ) {[Du, / Dt]y) with [Dui/Dt] v 
sampled homogeneously in space over the field (see Fig. [5}. 
Equation (fl2l) based on the first order approximation, although 
qualitatively correct, fails to predict quantitatively the measure- 
ments for dp > Arj. In Figure [3] we have plotted the data points 
from HB— . The agreement with our data is excellent up to 
dp Arj, while the Lagrangian model shows a less pronounced 
decrease (approximately by a factor of two) for larger diame- 
ters. We will analyze the possible origin of these difference in 
section 4. 

In Figure [4] we show the measurements of the acceleration 
flatness F(ai) = (df)/((a 2 )) 2 normalized by the fluid accelera- 
tion flatness F(Dui/Dt) versus particle size. As already noticed 
in—, and experimentally verified in— the flatness decreases 
with increasing the particle size. Here the different Lagrangian 
models lead only to small shifts in the flatness value, hence the 
picture remains the same as for the variance. Although, HB di- 
rect numerical simulations may suffers of lack convergence the 
case of flatness, the measurements are in qualitative agreement 
with our simulations. Finally we look at the correlation time 
of the particle acceleration. As done before ! °, we define the 
acceleration correlation time as the integral of the correlation 
function from time zero till its first change of sign: 



C a ,(r) dr ; C a ,(T) = 

a 



<a,(f + T)a,(Q) 
<(a,(0) 2 ) 



(13) 



where C ai (ro) = 0. We observe that, as soon as the particle 
size grows, the acceleration correlation time T a p deviates from 
the tracer value T a ,f(- 1.2 T n ), Fig. [5] This growth, is a re- 
sult of the Faxen averaging (in fact it is absent when averaging 
is not included, see the discussion in—) and comes from the 
fact that in finite-sized particles Stokes drag becomes negligi- 
ble, leading to dv/dt - [Du/Dt] v . We note however that there 
is also a significant difference between the basic Faxen Stokes 
drag model, for which the mechanism explained before is at 
work, and the model with non-Stokesian drag, the latter one 
producing more correlation. This feature is rather surprising 
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Figure 4: Single-component particle acceleration flatness, F(a,), normalized 
by the fluid tracer acceleration flatness, F(D,iij), as a function of the particle 
diameter in rj units. The dotted line represents the ratio of the flatness of the 
volume filtered eulerian field [Du/Dt] v to the fluid acceleration flatness. 



and has a different origin. One may think that including non- 
Stokesian drag, the effective response time of the particle, i.e., 
1/ T e// = (1 + CRe>,,)3vy8/r 2 is reduced, therefore the accelera- 
tion should be correlated on a shorter time scale ~ r e //. As 
we will see later in section [331 when non-Stokesian drag is ac- 
tive, the drag is never negligible and v - [u]s, hence for the 
acceleration dv/dt d[u] s /dt. It is clear that d[u] s jdt = 
dt [a]s + [u]s • <9 [u] s does not have the sub-grid (sub particle- 
size) correlations included in [Du/Dt] v = d t [u] v + [u • d u] v , 
which are correlated on shorter timescales. 
Any Lagrangian model equation seems however to underesti- 
mate the real T UtP resulting form the HB direct numerical sim- 
ulations. This fact has been also noticed, in a comparison of 
Faxen Lagrangian model with experimental data^. 

3.3. Velocity statistics 

Contrary to the acceleration's case the velocity particle statis- 
tics is deeply affected by the form of the particle dynamical 
equation. In Figure [6] we show as a function of the particle 
size the measurements of the deviation of normalized single- 
component particle velocity variance with respect to the ve- 
locity variance of the turbulent flow: ((u 2 ) - (v 2 ))/(u 2 ). The 
first model - based simply on Faxen terms and Stokes drag - 
predicts for this quantity a non-monotonous behavior, leading 
for the bigger particle sizes to a velocity variance even larger 
than the fluid one. This results is rather unphysical (as it is 
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Figure 5: Correlation time of acceleration T a p , in T n units, as a function of the 
normalized particle diameter d„/rj. 
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not possible for a particle to be on average more energetic then 
the flow by which it is driven and transported) and clearly it 
represents a limitation of the basic Faxen model with Stokes 
drag. This limitation is readily cured whenever an extra drag 
force is added. Among History and non-Stokesian drag it is 
definitely the latter one bringing the most significant changes. 
The non-Stokesian drag reduces the kinetic energy of the par- 
ticle as compard to the one of a fluid tracer, this energy de- 
creases monotonically with the particle size. For comparison, 
in Figure [6] we have also reported the prediction in the limit of 
small-particles, eq. ( fTTT ). As for the acceleration, this analytical 
prediction seems to be a good approximation to the measure- 
ments up to dp 4rj. We also report in Figure [6] the value of 
the variance of the Eulerian filtered field [a]s- One can note 
that by adding more drag to the basic Stokes term, the parti- 
cle velocity approaches the filtered fluid velocity variance, i.e. 
v - * [u]s- We note that the prediction of the Lagrangian mod- 
els keeping into account all the considered effects eqs. (H}-®, 
agree well with the HB data. Apart from the parabolic (~ dp) 
behavior for vanishing particle sizes no clear scaling of the nor- 
malized velocity variance can be detected form our measure- 
ment. In HB it was proposed that a scaling regime with slope 
of 2 / 3 would appear out of the Faxen dominated regime. This 
was attributed to the effect of the background turbulent flow, via 
the assumption (v 2 ) ~ ((dud p ) 2 } ~ dp where the Kolmogorov 
scaling relation for the second order eulerian structure function 
(((5m,.) 2 ) = ((u(x + r) - u(x)) • f 2 ) ~ r 1 ^ is implied. Here we 
would like to advance another explanation, the different scaling 
in that intermediate regime seems to be an effect of the drag 
term and in particular of the non-Stokesian drag (see again Fig- 
ure We note that non-Stokesian drag term included in our 
model equations based on SN parametrization accounts rather 
for the effect of a stationary wake behind the particle than for 
wake generated turbulent fluctuations. Given the good agree- 
ment of our data with HB measurements, our guess is that the 
background turbulent fluctuations plays only a minor role in de- 
termining the particle velocity statistics. 



Figure 6: Deviation of the particle velocity variance from the fluid value, as 
a function of the dimensionless particle diameter d p /rj. The behavior of the 
four different models adopted is reported. The dashed-dotted line represent 
the deviation from the fluid root-mean-square (r.m.s.) velocity that is expected 
to originate from Faxen corrections in the small particle limit (see eq (IT} ). 
The dashed line represents the deviation of the variance of the surface filtered 
eulerian field [u] s from the unfiltered velocity variance. 



The non-Stokesian term has also important consequences on 
the velocity correlation time. Such a time can be defined as the 
time integral of the correlation function: 



Jo 



C v , (t) dx ; C Vi .(t) 



<V,(f + T)Vf(f)) 

<(v,(0) 2 > 



(14) 



As we have already mentioned, the Stokes drag force alone is 
not effective in slowing down the velocity especially for large 
particles. This produces quasi-ballistic trajectories, governed 
by v(f) ^ C [D,u{t')] v df, that tend to have very large correla- 
tion time. The non-Stokesian (or wake) drag provides instead 
a way to reduce particle speed, v(f) ^ [u] s and its correlation 
time. This is evident from Figure [7] where the correlation func- 
tions for different terms for a large Stokesian and a large non- 
Stokesian particle are compared. 

On Figure [8] the trend of the velocity correlation time as a 
function of the particle size is shown. Non-Stokesian drag pro- 
duces a reduction of T p of more than 100% as compared to the 
purely Stokesian case. 

4. Discussion 

In previous sections we have shown that the agreement of the 
Lagrangian model, even in its complete form eqs. ([I]).!®,©, 
with the measurements from HB direct numerical simulations 
is rather satisfactory for the variance of the velocity as a func- 
tion of the particle size. However, we observe systematic de- 
viations when acceleration is concerned. Here we would like 
to discuss the causes of these discrepancies more in detail. We 



6 



CM t 
> 

A 



V 

\\ 




\\ Faxen Model with Stokes drag 

\ \ H 

\ \ 


:: \ v ^ 

\ \ 
V \ 

\ \ 
\ \ 




\ \ 

+ \ \ 













0.9 




0.8 




0.7 


CM 1 
> 


0.6 






P 


0.5 






> 


0.4 






V 






0.3 




0.2 




0.1 








|\ + non-Stokesian drag 

3 \ 

|.\ 


...0 D t U "')lv dt ' 


11 \ 

.11 \ 

±\ \ 
±\ \ 
±\ \ ■ 
; \ \ 

U j 

■ | \ \ 

I I : \ 

II \ 

_ .±. . \ \ 

I \ v 
± \ s 
± \ \ 

\ \ ; .% 

± \ N 
i \ ^ 







15 
i/T, 



Figure 7: CoiTelation function of particle velocity v, of the filtered fluid velocity 
along the particle trajectory [u]j , and of the time integral of fluid acceleration 
L[DtM(f)\yM for particle size d p = 32;/. Faxen model with Stokes drag 
(top), with non-Stokesian drag (bottom). 
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Figure 8: Integral correlation time of particle velocity T„ as a function of the 
diameter d p . T„ is made dimensionless by the integral velocity correlation time 
of a Lagrangian tracer Tl, while the diameter is normalized by the dissipative 
scale Tj. 



can advance the following hypothesis: 

i) The differences originate from a limitation of the model 
which takes into account volume and surface averages only in 
the approximate form of a Gaussian convolution. Although this 
approximation is well tuned for the first order Faxen correc- 
tion, it might be less accurate for larger particles, when roughly 
dp > 4rj. 

ii) The model neglect the effect of particle interaction with its 
own wake. This effect, while negligible in a flow with a large 
mean component, might become relevant in the isotropic flow 
conditions considered here (and in HB work) where a particle 
can cross a region previously perturbed by its own wake. 

iii) The observed discrepancies may come from differences in 
the simulated turbulent flows - in fact differences in the forc- 
ing at such small Reynolds number Re^ ~ 32 may have conse- 
quences even on the small scale statistics. The point i) can be 
investigated carefully. In order to see if there is any bias induced 
by Gaussian averaging as compared to the mean over a sphere, 
we have averaged the a field Du(x, t)/Dt over a large number 
of spheres of diameters d p , uniformly distributed over random 
locations in space. This procedure is repeated over 250 Eule- 
rian Du/Dt snapshots, equally spaced in time over an interval 
of 20 Te- Although this method is not precise for small sphere 
diameters when only few points of the discretized Du/Dt field 
enter into the average, for large diameters the average converges 
rapidly to the correct (continuum-limit) values. In figure [9] we 
report the results of these measurements for the variance and 
the flatness of [Du/Z)f]y (sphere), as a function of d p and we 
compare it with the Gaussian averages [Du/Dt]v (Gauss), and 
with the measurements on the Lagrangian models and DNS 
data. The result shows that the sphere average is sensibly dif- 
ferent from the Gaussian convolution. For small particles the 
discretization bias fails to match the analytical prediction <TT~TT > . 
which is instead well captured by the Gaussian distribution. For 
larger particle sizes the sphere average shows a stronger de- 
crease of the variance as compared to the Gaussian case. Re- 
markably in this region the average over spheres agrees with 
the HB results. A similar scenario appears for the flatness fig. 
[9] (bottom). This finding is important at least for two reasons. 
First it shows that the Faxen correction to the fluid acceleration 
plays a central role in the particle acceleration statistics not only 
for a particle in the upper dissipative range d p < \r\ but also for 
a particle with size is in the inertial range d p < 32rj. Second, it 
shows that, although the Gaussian convolution approximation is 
an efficient method to solve particle dynamics, it has limitations 
which becomes important approximately when d p > 4rj. The 
above observations have also an impact on the time-statistics: 
since a reduction of the variance of the filtered acceleration is 
associated to a slower fluctuations in time, we expect the corre- 
lation time of acceleration to increase for finite-sized particles. 
Hence the picture in fig [5] would change, solving the observed 
mismatch towards HB simulations and experiments^. How- 
ever, this latter point may find a confirmations only in further 
studies. 

Evaluating the impact of point ii) is unfortunately not possible 
in the framework of the present model. One would need to in- 
troduce a coupling (so called two-way coupling) between the 
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Figure 9: Comparison between acceleration variance (top panel) and flatness 
(bottom panel) We show data coming from Gaussian convolutions, same as in 
fig. [3] of the fluid acceleration field Du/Dt (dotted line), the quantity has been 
computed run-time therefore statistical errors are in this case of the order of 
the line thickness. Averages over spherical volumes (dashed line) are affected 
by a larger statistical uncertainty (yellow shaded region) which comes from 
the differences in the measurements between the three cartesian components. 
Measurements from the Lagrangian model based on eq ^ and eq. Q,j2).0 
are reported, together with HB data. 



particle and the fluid enforcing conservation of the total mo- 
mentum. 

Finally on point iii), for completeness we shall note that there 
are some differences between the turbulent flows simulated by 
Homann & Bee and the one used here. While we adopt a 
forcing which keeps constant the energy on the largest Fourier 
modes (amplitude-driving), HB also uses a random forcing of 
the phase in Fourier space (phase-driving). While these dif- 
ferences have no effect in fully developed turbulent flows, at 
the small Reynolds number considered here, Re,\ 32, they 
might have an impact. We see indeed that while we find 



-3/2 y l/2 
-3/2..1/2 



(af)s 
(af)s 

more intermittent at small-scales. This prevents us, for instance, 
from a direct comparison on the shape of the probability density 



= 1.1, F(cii) — 5.7, Homann & Bee report 
= 1.3 and F{a,) = 8.4, hence a flow slightly 



functions of acceleration and velocity. 



5. Conclusions 

In this study we have focused on the statistical properties 
of acceleration and velocity of finite-sized neutrally-buoyant 
particles driven in a moderately turbulent homogeneous and 
isotropic flow. We have adopted a Lagrangian model parti- 
cle equation which keeps into account inertia effect, size ef- 
fects, and the drag forces resulting both from a Stokes flow 
around the particle and form an asymmetric trailing wake state 
via Schiller-Naumann modeling. We have studied the contri- 
bution of these forces separately, in particular the drag force 
has been divided into three components: Stokes, History, and 
non-Stokesian force. 

We find that the drag forces have minor effect on the statis- 
tical properties of acceleration. Acceleration statistics seems 
to be dominated by inertia effect and by the Faxen corrections, 
whose influences extends also over particle with size in the in- 
ertial turbulent range. 

On the opposite drag forces have important effects on the 
time integral of the acceleration, that is to say on the veloc- 
ity statistics. This is particularly evident when considering the 
trend of the second order statistical moment (the variance) as a 
function of the particle size. For the case of neutrally buoyant 
particle analyzed here, the variance of the particle velocity from 
the different Lagrangian models start to separate at d p > 877, 
corresponding to a particle Reynolds number Re p ~ 0{lQ). 
Above this threshold History and non-Stokesian drag have a 
dominant role. This lead us to propose that the trend observed 
for the particle velocity variance as a function of its size does 
not originate from inertial-scale properties of the background 
turbulent flow but arise from the non-Stokesian component of 
the drag produced by the wake behind the particle. 

The effects detected in the velocity statistics are relevant for 
studies of particle dispersion in turbulence. For instance a sim- 
ple finite-size particle model based on Stokes drag or history 
force only, would overestimate the average particle dispersion 
from a fixed source in space, and similarly for pair separa- 
tion. Therefore, this study suggests that in order to validate La- 
grangian models one should look not only into the small scale 
acceleration statistics - as done up to now in many studies - but 
also into velocity and possibly dispersion and pair separation 
statistics. 

On the numerical side, we have shown that the Gaussian con- 
volution approximation, despite its computational efficiency is 
not accurate when particles much larger than 4rj are involved. 
On the other hand spherical averages in real space, as the one 
performed in our test, would be computationally very expen- 
sive and not enough accurate for small size particles (in partic- 
ular it is not possible to capture the laplacian correction ( fTUb 
when only few grid points are used). Future numerical studies 
should find a trade off between efficient computations and accu- 
racy. A possible way, which however needs careful scrutiny and 
tuning, is to consider the implementation of convolution kernel 
functions with sharper boundaries. Furthermore, the effect of 
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two-way coupling in this type of homogeneous and isotropic 
turbulent flow deserve to be studied. 
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